Flat band electrons and interactions in rhombohedral graphene multilayers 
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Multilayer graphene systems with a rhombohedral stacking order harbor nearly flat bands in their 
single-particle spectrum. We propose ansatz states to describe the edge-localized states of flat band 
electrons. The absence of kinetic dispersion near the fermi level leaves the interaction as a dominate 
mechanism to govern the low energy physics of a low density electron system. We build up an 
effective lattice model in two interacting low-energy bands, where the full terms of the Coulomb 
interaction, including those long-range and off-diagonal parts, have been considered. The interaction 
matrix coefficients in the many-body Hamiltonian model are directly calculated using orthonormal 
Wannier basis. We then present a flat-band projection to yield an interaction-only lattice model 
for flat band electrons. We find that this limited model might energetically favor a ferromagnetic 
quantum crystal. 

PACS numbers; 71.10.Fd, 71.27.-|-a, 73.21.Ac 



I. INTRODUCTION 

Graphene based structures have drawn numerous at- 
tentions due to their unique electronic propertieSfii^ The 
rapid technique development enable people to engineer 
the graphene nanostructures in special designs, yield- 
ing rich band structure features. In recent years, great 
theoretical--" and experimental^^"— interests have been 
focused on the graphene multilayer systems. Different 
from the graphene monolayer, the band structure of the 
multilayer graphene system depends on its stacking or- 
der, i.e., the way to stacking the graphene sheets. Re- 
cently, the rhombohedral stacking multilayer graphene 
has drawn intensive research interests due to its intrigu- 
ing band dispersion. It has two subbands near the neutral 
system Fermi level, one conduction band and one valence 
band with |e| ~ dispersion touching at e = 0, where 
N is the layer numberi^ The rather flat energy bands 
near e = make the rhombohedral stacking multilayer 
graphene susceptible to the interaction. Thus, the sys- 
tem is instable towards quantum correlated phases, such 
as superconductors or ferromagnets.^'^'^'* 

Some recent experiments^^'^^ in rhombohedral stack- 
ing graphene trilayer have shown the hints of a gapped 
ground state, which is in sharp contrast with the gapless 
semiconducting ground state suggested in noninteract- 
ing picture. Several symmetry-breaking correlated states 
have been proposed as the candidates of the gapped 
ground state, such as layered antiferromagnetic state, 
quantum anomalous hall state, quantum spin hall state, 
and quantum valley hall stateJ^!^ However, the theo- 
retical predictions strongly depend on the model and pa- 
rameters they chose. The detail properties of the ground 
state are still under debate. 

Flat band electrons of the rhombohedral stacking 
graphene system are of particular interest, since it is be- 
lieved that the correlated ground state results from the 
interplay between the electron-electron interaction and 
the peculiar flat energy bands near the Fermi level. For 



a low density system the dispersion-less flat bands leave 
the Coulomb interaction predominantly rule the low en- 
ergy physics. This calls for a comprehensive evaluation 
to the effects from all interaction terms, including those 
long-range density-density repulsion terms and leading 
off-diagonal terms, such as the direct spin exchange. The 
absence of the intra-band screening in a flat band sug- 
gests that these nonlocal interactions would be relevant. 
Studies have shown that these nonlocal interactions can 
lead to exotic correlated phases, such as quantum crystal 
and quantum liquids^S'^ In this paper, we theoretically 
investigate the flat band electrons and their interaction in 
the rhombohedral stacking graphene multilayer system. 
We establish a set of many-body Hamiltonian models, 
which allow to appropriately include the effects from non- 
local interaction in addition to the Hubbard onsite term. 
Corresponding to the unique non-interacting band struc- 
ture, a single-particle basis of Wannier functions is first 
constructed. We then use our basis to directly compute 
the matrix elements of a unscreened Coulomb interaction 
in two low-energy bands. A projection protocol has been 
presented to approach an approximate interaction-only 
lattice model in the flat-band limit, which are highly non- 
trivial, incorporating two bands, long-range interactions, 
and spins. We argue that, at low densities, the long-range 
part of the interaction in this limit model might support 
ferromagnetic quantum crystals. 

Our interaction model extends beyond the mean-fieldi^ 
and renormalization group ^^i^° studies, where a screened 
interaction with either the onsite Hubbard term or short- 
ranged interaction term is considered. Our study is also 
different from those with ab initio calculations^! and 
Hartree-Fock approximations which rely on certain lo- 
cal approximation to treat the nonlocal interaction and 
spin exchange terms. Alternatively, the Wannier basis 
allows us directly calculate these nonlocal terms. 

This paper is organized as follows. In Section [H] 
we consider the band structure that arises from the 
non-interacting tight-binding model of rhombohedral 
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FIG. 1: (Color online) Left: Schematic top-view of multiple 
layers of graphene sheets in the rhombohedral stacking or- 
der. Lines in solid, dot and dash types represent the in-plane 
carbon-carbon bonds at three neighboring layers, counted 
from top to bottom in the z direction. The shaded area cor- 
responds to a single unit cell. Right: Schematic side-view of 
a unit cell in a triangular prism shape. The solid and open 
circles stand for the atomic sites of the sublattice A and B, 
respectively, ji are three corner axes of the prism. 



stacking graphene systems. An ansatz wave function 
has been proposed to describe two flat bands. In 
Section |IlT] we construct localized single-particle basis 
states, orthonormal Wannier functions, from carbon tt^ 
orbitals in graphene honeycomb lattices. Section HVl uses 
the Wannier functions to explicitly compute Coulomb 
interaction matrix elements for two low-energy bands. 
Section |V] defines a projection scheme that limits the 
total many-body model to the flat-band portion of the 
single-particle spectrum and discuss the possible low 
energy physics of this interaction-only lattice model. 
Section I VI I summarizes and looks forward to more 
accurate studies of the models constructed here. 



II. FLAT BANDS IN RHOMBOHEDRAL 
STACKING GRAPHENE SHEETS 

We consider interacting electrons hopping among car- 
bon sites of rhombohedral graphene layers. In the left 
panel of Fig. [1] we schematically show the lattice of this 
stacking system. Two neighboring graphene layers have 
a relative in-plane shift along the carbon-carbon bond di- 
rection with the shift distance equal to the bond length 
Rq ~ I.42A. After three successive shifts, the forth 
layer recovers the same lattice as the first layer. Thus, 
Lz, the total number of stacking layers is a multiple of 
three. The layer separation is similar as the graphite with 
i?_L ~ 3.35A. As shown in the right panel of Fig. [1] the 
primitive unit cell is in the shape of a triangular prism 
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FIG. 2: (Color online) The dot-dashed lines indicate the en- 
ergy eigenvalues of Eq. ([T]) versus wavevector for the rhom- 
bohedral graphene trilayer. The solid line shows the approxi- 
mate expression for the energy, Eq. ([31 . Two flat bands form 
near the valley points K and K' . In the large limit, the 
bands flatten. 



with the total number of atom sites M = 2Lz ■ Each layer 
of the unit cell contains two sublattice sites of A and B 
with perpendicular bonds to their counterpart sublattice 
site at the neighboring layers. The array of unit cells 
forms a two-dimensional Bravais triangular lattice with 
the lattice length Rc = VSRo- 

In a simple non-interacting picture, the minimum 
single-particle tight-binding Hamiltonian is given ass^ 
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h.c), 



(1) 



where the sum is along carbon-carbon bonds and the hop- 
ping integrals are taken iy — 3.16 eV and t± — 0.39 eV 
for the intralayer and interlayer hopping, respectively^^ 
The second-quantized operator cjj creates a fermion at a 
site n. Labels n and m indicate lattice sites, in contrast 
to labels for unit cells, i,j,k,l, used in the following. 

Two bands near the Fermi level flatten around the cor- 
ners (K and K' valley points) of the Brillouin zone (BZ). 
An example band structure for a trilayer system is shown 
in Fig. 2. Crossing the Fermi level, the conduction band 
(upper band, u) and valence band (lower band, d) are 
nearly degenerate with in-plane wavevectors q (relative 
to the valley points) in a region |q| < and form flat 
bands. For larger number of stacking layers these bands 
can flatten considerably. 

To model the two flat bands and examine the band 
width, we construct analytical ansatz states in the lin- 
ear combination of atom orbital basis as {(f)A,(f>B)'^ 
with (l)A/B = (0A/B,i,=i(q), ••■,0A/s,j,=L,(q)), where 
the sites of sublattice A (B) on the bottom (top) layer 
have direct link to the neighboring layer. The indices iz 
marked from 1 to represent the graphene layers from 
the topmost one to the bottom as shown in Fig. 1. 

For a wave function to be exact for i? = 0, the math- 
ematical necessary condition requires the wave function 
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components between the neighboring layers to meet a 
certain relationship of 



p(q) = 



*A,i,+i(q) 



(q) 



'>B,t, (q) 
2 



_A[e-^i^Ro + 2cos(^(Z^i?o)e^^-^°/']. (2) 



Note that at the valley points of the K and K', we 
have |p(q)| = 0. The wave function is completely 
localized at two edge layers with the top layer occu- 
pied solely by the lattice A and the bottom layer oc- 
cupied solely by the lattice B. When the momen- 
tum is shifted away from the valley points, the wave 
function extends to the inner layers from the two edge 
layers. The ansatz wave functions in the vicinity of 
the valley points have the analytical form of <l'±(q) = 
{(j)A,±(pBV with (j)A = (l,_p(q), ...,_p(q)^'-~^) and 0b = 
{{p*{q))^'~^, ...,p*{q),l). In the general case with 
|p(q| ^ 1, this ansatz state associates with a non-even 
occupation of the two sublattice sites on edge graphene 
layers. 

Considering semi-infinite stacking layers of sublattices 
A (edge at the top surface) and B (edge at the bot- 
tom surface), the convergence of the wave function re- 
quires |p(q)| < 1. This determines the valid range of 
flat-band ansatz wave function with a radius g'^/IKj « 
{t±/t\\){V3/2n) in the limit of tj_/t\\ < 1. Here we see 
that an enhancement of the interlayer hopping leads to a 
larger flat-band sector. 

With the above ansatz states, the energy dispersion of 
bands T — u, d in the flat-band region can be computed 
explicitly: 
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l*±(q)^i^o(q)$±(q) 
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l$±(q)P 

|Re[p(q)^=]|(l 



b(q)l^ 



with 



Ha{q)=tj 



i-b(q)P^^ 

Q(q) \ 
Q^q) ) 



(3) 



indicating that the band dispersion plays a small role 
with the stacking number increasing. Such a vanishing 
bandwidth leaves the interaction as the dominant term 
in the full many-body Hamiltonian of electrons. 

For a dilute system with partially filled lattices, the 
lower-energy physics of the electron system is mainly de- 
termined by the single-particle basis states within the 
flat-band sectors near the Fermi level. Thus, we project 
the Hamiltonian into the basis of flat-band states in the 
approximation that Hq adds an overall constant energy 
shift to the spectrum. Our Hamiltonian model becomes: 



Hiatal = Y. ^r(q)c|j^rVr + -^^- 

q6BZ,o-,r 

constant + Vl^HyVp^ , 



(5) 



where the first equality is written in terms of the creation 



(annihilation) operator 



for a Bloch state at 



the wavevector q and band F, which is related to the 
operator for a single-particle basis state in the real space 
by the Fourier transform: 
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N ^ 

qGBZ 



'-qo-r- 



(6) 



Here Rj is the lattice vector of the j-th unit cell, N 
denotes the number of unit cells in the system and q- 
space mesh in the BZ, and a e {t, i} labels spin. "Ppg 
denotes a projection into flat bands such that the many- 
body eigenstates are constructed from Bloch states in the 
flat-band sectors |q| < . 

To explore possible many-body ground states in the 
rhombohedral stacking graphene system, we need con- 
struct an accurate form for Eq. ([5]) in the flat-band basis. 
The absence of dispersion excludes intra-band screening 
as in ordinary Fermi liquids.— Thus, many-body eigen- 
states are determined entirely by the interplay between 
various terms in the interaction. It is therefore crucial to 
accurately determine the interaction terms in Eq. ([SJ as 
prescribed by our choice of single-particle basis. In the 
next section, we describe how to construct orthonormal 
Wannier functions to serve as single-particle basis states. 



Q(q) 



/-p*(q) .. 
1 -p*(q) .. 



V 



1 -P*(q) 



As shown in Fig. 2, the analytical dispersion Eq. ([3]) 
agrees with those calculated directly from the tight-bind 
Hamiltonian in the vicinity of valley points, indicating 
the ansatz wave function as an effective approximation 
to flat-band states. 

With the equation ([3]) , we can estimate the bandwidth 
of the two nearly flat bands using the energy value at 
the flat-band boundary q"^. In the large Lz limit, the 
bandwidth for states in the flat-band sector vanishes as: 

A , , t \ 



|i?(q- 



(4) 



III. SINGLE-PARTICLE BASIS STATES: 
LOW-ENERGY BAND WANNIER FUNCTIONS 

In this section we superpose overlapping carbon tt^ or- 
bitals to form orthogonal Wannier functions. The Wan- 
nier functions will then be used to accurately determine 
interaction matrix elements in later sections. 

In an isolated band the Wannier functions are given by 



I^,(r) =I^o(r-R,) = 



(7) 



where momenta q sum over N discrete values in- 
side the entire BZ. The Bloch functions are ^q(r) = 



X/m=l C'mqXiTiq(r). 
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summation over all local atomic orbitals <^(r), i.e. 
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FIG. 3: (Color online) Wannier functions of trilayer graphene 
sheets, (a) u-band case: Two three-dimensional plots on the 
left represent the distribution of the Wannier functions in the 
xy plane at the 2-positions right above the top layer and right 
below the bottom layer, respectively. The twin peaks locate 
at two sublattice sites of the edge layer in the original unit 
cell. The cartoon on the right plots the distribution of the 
Wannier function in the z direction along three corner axes 
of the original unit cell, (b) The same plots as (a) but for 
d-band case, (c) u-band Wannier function with t± = ty. 



To make the contact with first principles calculations^^ 
we form Bloch functions from carbon tt^ orbitals, 0(r) = 



(i/ViV)Ec 



The 



basis states become 



Xmq 

r\- -rmj), where r^j = -I- T 
the location of the m-th atom in the j-th unit cell. 



(r) = 
is 



The coefficients C„iq and energy eigenvalues E{q) are 



obtained from diagonalization of the secular equation: 



(8) 



Woir) = Nf 



M 

EE' 



(9) 



with weights a„jj — Cmqe'^i'^J and the normaliza- 
tion constant Nf. A denser sampling in momentum space 
(i.e., larger N) yields more accurate Wannier functions. 
In practice, we find that the Wannier function has al- 
ready converged when taking N = 1261 for — 3. 

We can extend our construction of the Wannier func- 
tions to include both the upper and lower bands. The 
Wannier functions of these two low-energy bands in a 
trilayer system are shown in Fig. [3l We note that these 
two Wannier functions mainly localize at the original unit 
cell with the reflection symmetry (antisymmetry) along 
a center line (-s/S, 1,0) for the upper (lower) band, de- 
caying rapidly within several cell lengths. In the plots of 
Wannier functions as a function of z-positions as shown 
in Fig 3. (a) and (b), Wannier functions mostly distribute 
in a narrow region around each layer with the node on 
the layer. This is due to the property of underlying tt^ 
orbitals. We note that there exists a large portion of the 
Wannier function around the middle layer, indicating the 
contribution from those extensive states with momenta 
outside the flat-band region. Under the given hopping 
parameters of t_L/t|| ^ 0.1, two sublattices near evenly 
occupy each layer. 

To explore the effect of the ffat-band on the Wannier 
functions, we study the case with t±_ = t\\. As shown 
in Fig. 3(c), the non-balanced occupation between sub- 
lattices A and B at two edge layers magnifies as the 
ffat-band region expands, consistent with the property 
of the ansatz ffat-band state in the previous section. 
Meanwhile, the relative portion of the extensive Wan- 
nier function around the middle layer also reduces as ex- 
pected. The ffat-band induced asymmetric occupation of 
two sublattices in the edge layers may justify the gapped 
symmetry-breaking states proposed by earlier theoretical 
studies.—!^ 



IV. COULOMB INTERACTION MODEL 



where the matrix H follows from the tight- 
binding Hamiltonian Hq with elements H(q)mn — 
J (irXmq(r)-ffoXnq(r). The overlap matrix O are taken 
as the the identity matrix in the tight-binding approx- 
imation. The eigenvectors Cq = {Ciq, Ci/q}^ yield 
the coefficients used in the definition of the Wannier 
functions. 

To construct orthonormal Wannier function, a specific 
set of single-particle basis states are chosen by enforcing 
Cmq at the edge atomic sites m = 1 and m = M con- 
jugate. The resulting Wannier functions Wj{r) are real 
and localized at Rj with the certain symmetry between 
top and bottom portions in the stacking direction z. 

We can write the Wannier function at the origin as a 



For a dilute system where the chemical potential lies 
between the two nearly fiat bands, the Coulomb interac- 
tion can in principle favor occupancy of both bands or a 
single band. As a ffrst approximation, we assume that 
the valence band is inert and only the conduction band, 
u, is active. 

An unscreened Coulomb interaction in a single band 
has a second-quantized many-body form of 



z,j,fc,/;(T,cr 



/C, 
ja ha 



(10) 



where the operators c]^ (c,^) create (annihilate) a 
fermion with spin cr in a Wannier state centered at the 
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TABLE I; Fitting parameters for the Gaussian approximation 
to the TT^ orbital with ^ — 1.72. 



s 


1 


2 


3 


Is 
Ps 


0.15591627 
2.9412494 


0.60768372 
0.6834831 


0.39195739 
0.2222899 



TABLE IL Matrix elements for one-band (u band) case of the 
1/2=3 system with unit cell separations of up to 3i?c. 



Vo=3.587e-l 


1 R-i — R-j / Rc 

Jij 
Vj 


1 

2.136e-3 
2.007e-l 


2.232e-3 
1.462e-l 


2 

9.703e-4 
1.319e-l 


V7 
5.273e-4 
1.080e-l 


3 

6.075e-5 
9.751e-2 



TABLE III: Same as the Table lU but for the two-band case. 



Vq — i.4yi3e-i 




Vo"=3.587e-l 


T/ 1 OOq 1 

— z.iyze-i 




^0=4.9326-2 


j^j— y.c5D4e-z 








jJVi — ±Vj j/ltc 


-1 

i 


2 


3 




o /1 1 o 
0.4iye-Z 


5.527e-2 


4.079e-2 




2.007e-l 


1.319e-l 


9.751e-2 




i.oU4e-i 


8.562e-2 


6.313e-2 


■J I] 


8.877e-4 


3.965e-4 


2.443e-5 




2.136e-3 


9.703e-4 


6.075e-5 




2.726e-4 


2.052e-4 


2.058e-5 




7.467e-4 


7.920e-5 


2.994e-5 




6.885e-4 


3.100e-4 


3.853e-5 



i-th unit cell. The matrix elements V are determined by 
the Wannier basis given in last section. We can rewrite 
the above many-body Coulomb interaction in a sugges- 
tive form: 

i i<j i<j 

+ I E (11) 

{ij}(2{fe,;};<T,<T' 

Here, the single-component and total density operators 
are rii^ — cl^c^^ and = rii-^ + Ui^, respectively. The 

spin operators Si — (1/2) X^o-o-' '^la^aa' '^i^' defined 
in terms of the Pauli matrices ct. 

Eq. ([Tl]) keeps all terms in the full Coulomb interac- 
tion. The first term is the ordinary onsite Hubbard term 
used in some mean-field studies of multilayer graphene 
systems.—tiS- The second term captures the diagonal por- 
tion of the Coulomb interaction at long range, which 
favors certain charge order, such as a two-dimensional 
Wigner crystal. The absence of a dispersion in a low 
density system implies that this term can be relevant and 
must be kept in accurate models. The third term, the di- 
rect exchange term, favors ferromagnetism for Jij > 0. 
The last term represents remaining off-diagonal terms 
due to the Coulomb interaction, which are much small 
compared to the first three terms for a single band ac- 
cording to our direct calculation. 

The matrix elements in Eq. (jlip can be computed ex- 
plicitly using the Wannier basis in the u band, as shown 
in the appendix, Eqs. (fT6|). To calculate these inte- 
gral equations, we have approximated the exponential 
part of the tt^ orbital as a linear combination of three 
Gaussian functions: (/.(r) « 7^(128/3f/7r3)i/4ze-^=''', 
where the parameters 7s and /3s are obtained from the 
ST0-3G package."*^ Data for fitting the tt^ orbital with 

= 1.72 are listed in Table HI For numerical results shown 
here and in the following tables, we use the Bohr radius, 
ao = O.53A, as the unit length and the Coulomb energy 
e^/Aireao 27.2 eV in vacuum) as the unit of energy. 

Table |lT] lists the coefficients computed for an example 



system with = 3. As we see, all coefficients are pos- 
itive and can be sorted by Vq > Vij > Jij > 0. These 
coefficients suggest that a partially filled single band sup- 
ports the formation of ferromagnetic crystals. 

However, the large Coulomb interaction may cause 
mixing between the u and d bands. We need consider 
a more comprehensive two-band interaction model with 
Wannier functions in both the u and d bands. The inter- 
action Hamiltonian is dominated by the following terms: 

*,r i,a,r^r' 

t \r^r' J 
+ E (^^J "^mjT - JfjS^r ■ Sjr) 

+ E E (^u"«r?ijT' - J'tj^tr ■ SjT') 

+ E E E(^y^«T£7cJrvSr^' 

+ ^JckeJrV'VySrJ- (12) 

We have checked, by direct calculations, that other terms 
involving three or four centers are much smaller than 
terms kept here. In Eq. [T^ we see the Hubbard and 
ferromagnetic terms as in the one-band analysis. Besides, 
we have the non-trivial band exchange as the last term. 
The integral equations for all coefficients in Eq. [TH are 
listed in the Appendix. 

Table IHII shows numerically computed coefficients for 
the two-band model Eq. [12] in the example trilayer sys- 
tem. Rows 1-3 exhibit several leading terms of the diago- 
nal components of Coulomb interaction, which primarily 
determine the charge degrees of freedom. Rows 4-6 gov- 
ern the spin degrees of freedom. The positive elements 
support ferromagnetism. The last two rows give rise to 
band exchange effects. 

Our results support the in-plane ferromagnetic order 
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suggested by several theoretical models.J^"— The coeffi- 
cients of the onsite Coulomb repulsion have values of 2 — 5 
eV with an estimated effective dielectric constant e = 2 in 
graphene systems, consistent with the parameter range 
in a mean-field analysig^*^ for the experimentally observed 
energy gap4ii^ We also note that the long-range inter- 
action terms of up to the fifth nearest neighbors (rows 
1-3) have a magnitude comparable to the onsite terms, 
indicating the effective interaction range could be much 
longer than the usual screened interaction treatments 
with up to nearest or next-to-nearest neighbors. Based 
on the energetic argument these long-range terms are rel- 
evant and should be included to discuss the possible low 
energy states of a dilute system. 



V. FLAT-BAND PROJECTION 

In this section we construct a set of operators that al- 
low flat-band projection of the many-body Hamiltonian 
model constructed in the previous sections. We then dis- 
cuss the possible low energy physics under certain condi- 
tions. 

To enforce flat-band projection we limit all g-space 
sums to the flat-band region (FBR) |q| < q'^. We first 
consider a state operator in a single band that limits itself 
to the FBR: 

^l-^E E ')cL. (13) 

/ qGFBR 

This operator creates states centered around the unit cell 
at Rj while can overlap considerably with neighbors, in- 
dicating that the projection into a flat band delocalizes 
basis states. In the limit that the flat band encompasses 
the entire Brillouin zone, we have b^^ — )> c|^. 

We can then rewrite our model in terms of projected 
density and spin operators. The single-component and 
total projected density operators are pia = bl^b-^ and 
Pi = pi^ + pii, respectively. The projected spin operators 
are defined as: 

^, = — F F e'^''"'''^'^'cl,,a'C,,. (14) 

era' q.q'eFBR 

Note that these projected operators do not exhibit ordi- 
nary commutation relations because the underlying basis 
states are delocalized. 

The projected Hamiltonian can be rewritten entirely 
in terms of the above projected operators. Starting from 
an unprojected interaction model, we impose projection 
using the following replacements: c — ^ 6, n — > p, and S 

Considering the intrinsic energetic ordering as shown 
in the table IIIIl we rewrite the two-band interaction 



Hamiltonian in the projected space: 

+ E PiTPjT' - J^^f $iV ■ $jT') 

+ ^^Band-oxch, (15) 

where we have redefined the diagonal Coulomb terms: 

V\tf ^ Vlr y"^"'^'^^ ^ VL and V^=]' ^ other- 
— r,r' 

wise V ij — 0. We have also redefined the off-diagonal 
1 -rr#r' ^. — r=d.r'=M ^, , 

exchange terms: J i^j = Jlj, J a = J a, and 

p — p' p p' 

Jj^j = Jfj, otherwise J.^- — 0. The last term in 
Eq. (fTS)) corresponds to the last term in Eq. (|T2|) . 

Given a dilute system where the d band is inert and the 
band-exchange terms can be ignored, the first three terms 
in Eq. (jl5p will impose a rigid charge order analogy to the 
classical Wigner crystal. However, here the charges are 
significantly delocalized. A finite overlap among neigh- 
bors implies that the charges exist in a superposition of 
several different unit cells at once, indicating a quan- 
tum crystal. The forth term corresponds to an effective 
Heisenberg model which favors ferromagnetism between 
neighboring cell spins. Thus, the projected system fa- 
vors the ground state as ferromagnetic quantum crys- 
tals. Correspondingly, low energy spin excitations could 
be ferromagnetic magnons4^ 

VI. SUMMARY AND OUTLOOK 

We construct interaction lattice model for flat band 
electrons in rhombohedral stacking graphene layers. An 
ansatz wave function has been proposed to describe the 
properties of flat-band states emerging in the single- 
particle spectrum of the system. A single-particle basis 
of orthonormal Wannier functions was built from carbon 
TTz orbitals of the underlying graphene honeycomb lattice. 
We use this single-particle basis to explicitly compute the 
Coulomb matrix elements. The total model, Eqs. (|12l) . 
was then projected into the fiat bands, suggesting a fer- 
romagnetic quantum crystal ground state. Our flat-band 
model, Eq. P3t . sets the stage for more accurate study 
with a combination of variational studies and diagonal- 
ization to verify possible ground and excited states42i 

We want to stress the difference between the work pre- 
sented here and a previous mean-field study^''^ Our inter- 
action model includes a full consideration of the nonlo- 
cal interaction terms from two low-energy bands. The 
mean-field study^^ takes interaction contribution from 
all bands but only counts the onsite interaction term. 
Our model can be applicable in the limit case with the 
large stacking number and weak interaction. There the 
low-energy properties of the system are most relevant to 
two extremely fiat bands. In an otherwise case where 
the interaction is strong and the screening effect from 
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those dispersive higher-energy bands are not neghgible, The last term is used only in Eq. 
the mean-field treatment would be justified. 

The constructed model focuses on key physics of in- 
teracting flat bands but excludes a number of realistic 
effects. The distortion from long range hopping and the 
experimental conditions, such as defects and substrate 
disorder can destroy the flat-band approximation. Inter- 
band screening from nearby bands could also lead to cor- 
rections to the Coulomb interaction studied here. 
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VIII. APPENDIX 

The coefficients in Eqs. [TT] and fT2l are given by: 

Vo^ - / |^ror(r)W^or(r')P, 

4, - ^1 ^^w:Ar)W,r{r)W,rir')W;Ar'), 

4 = ^1 ^^W:^ir)W,,ir)WUr')W;Ar'), 

^^J = J y—^W:^ir)W.,{r)W,^ir')W*Ar'), 
= I ^^WUr)W,^ir)W,4r')W*Ar'), 

y—^^W:MWUr)W;^{r')WU^'). (16) 
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